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The work described in this report was performed by the Earth and 
Space Sciences Division of the Jet Propulsion Laboratory for the 
National Institutes of Health, Division of General Medical Sciences, 
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ABSTRACT 


An electro-mechanical scanning assembly has been designed and 
fabricated for the purpose of generating an ultrasound tomogram. A 
low cost modality has been demonstrated in which analog instrumenta- 
tion methods form a tomogram on photographic film. Successful tomo- 
gram reconstructions have been obtained on in-vitro test objects by 
using the attenuation of the first path ultrasound signal as it 
passes through the test object. Thus, the nearly half-century tomo- 
graphic methods of X-ray analysis have been verified as being useful 
for ultrasound imaging. 
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Summary 

An electro-mechanical scanning assembly has been designed and 
fabricated for the purpose of generating an ultrasound tomogram. A 
low cost modality has been demonstrated in which analog instrumenta- 
tion methods form a tomogram on photographic film. Successful tomo- 
gram reconstructions have been obtained on in-vitro test objects by 
using the attenuation of the first path ultrasound signal as it 
passes through the test object. Thus, the nearly half-century tomo- 
graphic methods of X-ray analysis have been verified as being useful 
for ultrasound imaging. 

Good geometric reconstruction of simple objects has been obtained 
over a field of 6 cm by 8 cm. With a system bandwidth of 1 MHz, cen- 
tered at 2.5 MHz, and unfocussed transducers measuring 5 mm in diameter, 
a resolution of approximately 3 mm has been demonstrated by using the 
coherent signal processing technique known as time delay spectrometry. 
Computer simulation has yielded results sensibly in agreement with the 
analog experimental data. Thus, the original objectives of this grant 
have been successfully met. 

Additional considerations have been revealed by the many physical 
measurements which were made during the course of the task. First, al- 
though quite acceptable images are obtained by using photographic film 
as the reconstruction medium, the inability of film to subtract light 
values creates an average background bias level on the tomogram. Dis- 
crete element memory matrices show promise of replacing film as a low 
cost reconstruction medium. Second, the need for an adaptive time 
window for accepting the earliest sound has been noted for those ob- 
jects which are both highly elongated in geometry and exhibit substan- 
tial deviation of acoustic velocity from the surrounding medium. 
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Background of Tomogram Technology 


It is a well known fact that it is possible to infer the structure 
of a multidimensional object from a set of its projections (see 
Appendices A and B). The mathematical basis for this reconstruction dates 
back to 1917 (1). The first practical use of this to medicine is due to 
Ziedses des Plantes (2), who formulated the basic principles of X-ray 
tomography. 

His work pioneered what is now called body-section radiography, and 
his original techniques proved quite successful in clinical application. 
The des Plantes process produces images by preserving the differential 
attenutation profile of a planar cross section. The images are built up 
on a photographic emulsion. An X-ray point source and a film cassette 
are articulated about the object in such a manner that all points on the 
tomographic plane are constrained to remain invariant on the film, while 
points not. on the tomographic plane are blurred by the scanning motion 
of film and X-ray source. As the scan progresses, the image of the tomo- 
graphic plane is continuously reinforced, while the other subject planes 
are blurred by the scanning motion. As simple as this concept might be, 
even considering the limitations of the usable range of film exposure, the 
results were considered clinically useful for more than four decades. 

What makes this technique work is the fact that X-rays principally 
react through differential absorption and seldom rebound to create sig- 
nificant diffraction of reverberation. Working against this technique is 
the weak interaction of X-radiation with soft body tissue and the addi- 
tive properties of film which build up both the out of focus as well as 
the focussed parts of the image. 

The advent of computerized axial tomography (CAT) has virtually 
eliminated these difficulties. Computer memory, rather than film, is 
the storage medium and signal values can thus be processed as desired. 

The weak interaction with soft tissue is not eliminated in CAT, since 
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this is a result of the physics of the process, but differential atten- 
uation can be extracted for more efficient image formation. 

Whereas, with des Plantes' original film method, the tomograph 
plane was essentially perpendicular to the line of centers between X- 
ray source and film, CAT creates a tomograph plane defined by the X-ray 
path. This change of plane has been regarded as an improvement in CAT 
since the image format includes the outer boundaries of the object being 
viewed. 

There is one serious drawback with present CAT relative to des 
Plantes' original method, and that is its cost. The type of equipment 
required for scanning and the nature of computer reconstruction has 
elevated the cost of such equipment well beyond what a small clime 
might afford. In addition, the use of X-rays still poses a hazard, 
no matter how slight, to the patient who is scanned. Finally, the 
physics of tissue/X-ray interaction still present a limitation to the 
subtleties of tissue pathology which can be visualized in the presence 
of photon noise. 

Ultrasound, unlike X-ray photons, interacts strongly with soft 
tissue; and the interaction is of a nature quite different than that of 
photons. X-ray photons do not sensibly speed up or slow down in passage 
through tissue. But, ultrasound velocity is itself an indicator of one 
aspect of tissue property. The differential absorption of ultrasound 
is due to factors which are different than those which cause absorption 
of X-ray photons. Thus, it can be argued that even if X-ray CAT had no 
safety or cost problems, it is worthwhile to pursue ultrasound tomography 
since it visualizes other tissue properties. Ultrasound can thus comple- 
ment X-ray images. 

It was recognition of these factors which led to the work discussed 
in this report. It was decided to investigate the use of ultrasound for 
the generation of a tomogram. TKe format of this tomogram is to be of 
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the type now obtained with CAT. An ultrasound source and receiver are 
positioned on opposite sides of the object to be scanned. A method is 
used that selects the earliest sound through the object, which has the 
greatest likelihood of being the direct sound analogous to the X-ray 
path. This sound is processed to produce a signal proportional to the 
true energy of signal passage, regardless of waveshape deformation due 
to differential frequency attenuation of the intervening tissue. The 
reconstructed tomogram is then converted to a photographic image for ob- 
servation and analysis. 

It was the intent of this task to find out just how inexpensively 
it might be possible to generate an ultrasound tomogram, recognizing that 
to do so will produce a sub-optimal image. The question was, are the re- 
sults of the lowest cost tomogram still clinically useful? Or, in collo- 
quial language, if one removes all the bells and whistles, is it possible 
to envision a' class of instrument that could find use in small clinics? 
The result of this present investigation is an affimative answer to 
these questions. 
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Definition of Terms 


X-Ray and Ultrasound Tomography 

In medical terminology, the word tomograph traditionally referred 
to an X-ray photograph of a selected plane in the body. It had the geo- 
metrical appearance of a tissue cross-section. The advent of ultra- 
sound imaging using pulse-echo techniques led to a format that was geo- 
metrically identical to the X-ray tomograph in the choice of coordinates, 
but was an image of ultrasound reflection properties and not those of 
electromagnetic radiation. By convention, this type of image came to be 
known as an ultrasound tomograph. 

Up to the present time, this duplication of terminology could 
cause no confusion because of the different modalities involved. How- 
ever, the recent introduction of image reconstruction methods capable 
of generating a sectional view from transil luminated projections may 
already be the source of some confusion, particularly in the case of 
ultrasound. 

There are, at present, three basic types of tomograph images, one 
using X-rays and two using ultrasound. The information contained within 
these images is complementary. With very few exceptions the images will 
show different information. One is not a replacement for the other. 

The differences among the three tomographs may be seen by inspecting 
the three images of the same section of the body. One of these will be an 
X-ray tomograph and will show the differential attenuation of tissue for 
X-rays passing through that tissue. The second will be a conventi^onal 
ultrasound pulse-echo tomograph and will outline the boundaries between 
tissue of different acoustic impedance by indicating the amount of sound 
reflected back from those boundaries. The third will be an ultrasound 
tomograph made by reconstructing the information obtained by the passage 
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of ultrasound completely through the section. This portrays the dif- 
ferential attenuation for sound due to passage through various types 
of tissue. 

With the exception of the external boundaries of the subject and 
a few dominant structural characteristics, these three types of tomo- 
graph probably will not look alike. That does not mean that any one 
is better than another. Rather, it means that the type of information 
contained within them is different. Each kind of tomograph can reveal 
a great deal of information to someone experienced in interpreting that 
particular type of image. What is potentially more significant, how- 
ever, is the synergism that may occur with two or three types of tomo- 
graph, each revealing its own peculair set of information. This could 
lead to a diagnostic capability not available from the use of any single 
type of tomograph. 

The class of tomograph which this report deals with is that ob- 
tained from a measure of the ultrasonic energy that passes completely 
through the tissue. This is an ultrasound transmission (or attenua- 
tion) tomograph and, in a general sense, is the acoustic analog of the 
X-ray tomograph. Because of the more complete control and processing 
that can take place with ultrasound signals, a great deal more relative 
information is obtained from this ultrasound signal than is obtained 
from X-rays alone. 
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Design of Analog Equipment 


As proposed, an attachment was developed for converting the existing 
JPL rectilinear scanner to a tomograph scanner. The rectilinear scanner 
uses two precision stepper motor drives mounted at right angles. Trans- 
mitting and receiving ultrasound transducers are mounted on a rigid yoke 
assembly and articulated in a raster, scan fashion about a test specimen. 
The transducers are positioned so as to be on opposite sides of the test - 
specimen and are mounted such that the major axis of the two transducers 
are aligned for transillumination of the specimen. 

The scanner provides two degrees of translational motion, but does 
not provide rotational motion about the test object. Because of the 
substantial mass of the drive mechanism, it was not considered feasible 
to convert the drive to provide tomographic reconstruction. A simpler 
solution, and one consonant with the requirements of this task, was to 
rotate the specimen and allow the scanner to translate the ultrasound 
probes past it in such a manner that the ultrasound signal always re- 
mained in one plane through the specimen regardless of position. The 
manner in which this was done was outlined in Figure 5 of the proposal 
and is reproduced here as Figure 1 of this report. 

A multiple-position rotary switch of robust construction was procured 
and used for support of the specimen as well as an angle pickoff. A 48- 
position indexing switch provides repeatable angle measurements every 7.5 
degrees around the whole of the test object. 

Electrical signals corresponding to the sine and cosine of the angle 
of test object rotation relative to scanner axis are obtained from two 
quadrature wipers on the switch assembly. Precision resistor ladders are 
tapped by the wipers to provide the equivalent of a discrete sine/cosine 
potentiometer. 
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Fig. 1. Symbolic Representation of the Combined Physical Translation of 
Ultrasound Crystals and Rotation of Test Object Which is Used to 
Generate a Tomographic Image. 
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The purpose of the electrical pickoff is to allow a multiplicity 
of angular views to be taken and registered with respect to the probe 
scanning geometry. The angle pickoff is processed to provide a modified 
C-mode scan on a display tube. In the visual presentation, the coordin- 
ates of the display device are fixed to those of the object being scanned. 
As the object is rotated for subsequent views, the horizontal and verti- 
cal drive to the display tube is inversely programmed so that the recon- 
structed object remains firmly placed on the face of the tube. The re- 
construction geometry is, thus, such that a fixed point on the tomographic 
section corresponds to a fixed point on the face of the display tube. 
Arithmetic reconstruction of the tomograph is performed by utilizing a 
photographic place to "add" the spot intensity at each place on the 
eventual tomograph. The basis for this is the time-exposure reciprocity 
known to exist over a limited range on all photographic emulsions. 

A simplified diagram of the electronic circuitry designed for the 
display tube drive is shown in Figure 2. In order to simplify the dia- 
gram, only that circuitry necessary to provide the modified C-mode scan 
is shown. There was a moderate amount of additional circuitry required 
to interface with JPL's specialized equipment, but that is not essential 
to an understanding of the operation of the tomograph scanner. 

The signal referred to in this figure as the horizontal input is 
a voltage proportional to the translational position of the scanner. 

Two additional voltages are added to this translational position signal. 

A one- kilohertz sawtooth of fixed amplitude is added for the purpose of 
creating a straight line on the face of the display tube. This straight 
line will create the trajectory of the geometric path of the ultrasound 
signal passing through the test object. The second signal added to the 
horizontal input is the tomographic center control. This allows the 
operator to readjust the center of the cathode ray display to be aligned 
with the center of the test object. 
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These three signals are sent in anti-phase relationship to the mul- 
tiple-position rotary switch. In this diagram the 48-position switch is 
symbolized by the circle, denoting the series-connected resistors that 
connect adjacent terminals. There are two wiper arms mounted at ninety 
degrees with respect to each other. These are symbolized by the arrows. 

In order to provide a straight line on the screen which is in quad- 
rature with the direction of mechanical scan, and thus aligned with the 
ultrasound energy, it is necessary to provide the low-pass and high-pass 
circuitry that follows the rotary switch pickoff. This circuitry directs 
the low-frequency voltages of translational motion to one axis and the 
high-frequency sawtooth voltage to the quadrature axis of the display 
tube. 
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Test Results 


The general results of the experimental program are perhaps best 
presented as resultant photographic images. Throughout this program it 
was decided to use Polaroid film of the type commonly used to photograph 
ultrasound images. It was realized that such film, with its very high 
contrast, is a poor choice for use in reconstructing tomographic images. 
Nevertheless, if the concept worked with such film, it would demonstrate 
an extremely low cost modality. 

In all cases, an open shutter exposure mode was used and the final 
image obtained as a superposition of all scans. The cathode ray display 
was carefully adjusted so that the major reticle marks which are seen in 
these images corresponded to one centimeter calibrations in the object 
space. The test objects were carefully positioned to lie in various ran- 
dom positions relative to the geometric center of the tomogram reconstruc- 
tion. This was done in order to assure that no inadvertent image improve- 
ment could occur if, for some reason, the geometric reconstruction depended 
upon the relative position of the object and scanner. 

Both 360 degree reconstructions and 180 degree reconstructions were 
made on a number of test objects. No difference could be noticed between 
these views, which indicates the reciprocity of scan direction. 

One of the first sets of test objects which were run is shown in 
Figures 3 and 4. This consists of a cross section of a circular hypo- 
demic case and a 4-40 metal screw. The hypodermic case is one centi- 
meter in diameter and approximately 1.5 millimeters in thickness. Ex- 
cept at grazing incidence, ultrasound passes through the lumen of the 
case. The metal screw is approximately 2 millimeters in thickness and 
represents a point scatterer smaller than the anticipated resolution of 
this system. Transmit and receive probes are 5 millimeter diameter, un- 
focussed discs and are 22 centimeters apart. The frequency range is 2 
to 3 megahertz. 
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Figure 3. Positive Image of Hypodermic Case 



Figure 4. Rescan Image of Hypodermic Case Using Inverted Polarity 
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Figure 3 is a positive image and Figure 4 is a rescan using inverted 
polarity to the cathode ray intensity modulation. As anticipated, a point 
scatterer appears as a first-order Fresnel ring with about a 6 millimeter 
diameter. In Figure 3 the additive background is the result of using film 
which can only add, never subtract, light values. The geometric place- 
ment and configuration of the reconstructed image corresponds within a 
millimeter of the known position of these test objects. 

Figures 5 and 6 are images of the hypodermic case alone. In Figure 
6 the video signal to the intensity modulation terminal was high-pass 
filtered to determine whether edges could be sharpened by such simple 
means. The results are poor. 

Figures 7 and 8 are two views of a low-contrast and high-contrast 
image, respectively, of a large diameter (1.3 cm) artery with associated 
soft tissue. This specimen is formalin fixed. The tomogram section is 
taken at right angles to the artery axis, just as the previous geometry 
of the hypodermic case, and the circularity of this section is in evi- 
dence. The shadowy outline around the periphery of the specimen is the 
outline of the plastic bag which held this sample in the water bath for 
tomogram scanning. 

Figures 9, 10, 11, and 12 are images obtained on a kidney specimen 
which was formalin fixed and contained in a plastic bag. Figure 9 is a 
high-contrast transmission shadowgraph image obtained by the JPL raster 
scanning method. This image is the ultrasound analog of an X-ray shadow- 
graph and the method of obtaining such an image has been explained else- 
where. Figure 10 is a tomograph reconstruction of the density distribu- 
tion in the plane lying 2.5 centimeters upward from the bottom of Figure 
9. The small variation of sound energy through such kidney specimens 
causes an extremely low contrast in the reconstructed image. Neverthe- 
less, a recognizable tomogram results. Figures 11 and 12 are tomograms 
corresponding to different levels through this object. The purpose in 
showing these results is to illustrate that in-vitro tomograms of tissue 
with small density variations do appear possible. 
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Figure 6. High-Pass Filtered Reconstructed Image of Test Object 
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Figure 7. Low-Contrast Image of a Large Artery 



Figure 8. High-Contrast Image of a Large Artery 
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Figure 9. High-Contrast Shadowgraph Image by Raster Scanning 
Method of a Kidney Specimen 



Figure 10. Tomograph Reconstruction of the Density Distribution 
of a Kidney Section 
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Figure 11. Tomogram at a Different Level of the Kidney Specimen 



Figure 12. Tomogram at a Still Different Level of the Kidney Specimen 
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A more substantial variation in density produces stronger images, as 
shown in Figure 13. This specimen is one which has been used by this ex- 
perimenter in a much earlier paper. It is a one-centimeter thick brain 
section containing two glioblastoma, which are distinguished under ultra- 
sound analysis by an unusually high attenuation. The shadowy outline of 
Figure 13 is due to the bag containing the specimen. The small differen- 
tial attenuation within the bag relative to the medium outside the bag 
is believed to be due to the interface attenuation between the bag and 
formalin. The edge of the brain section does not show in the tomogram 
reconstruction, possibly because of the masking effect of the shallow 
grazing angle, which the edge of the bag makes with those ultrasound rays 
parallel to the edge of the brain. However, the glioblastomas are well 
outlined and appear as the darker objects in the lower right quadrant of 
the image. 
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Computer Simulation 


Simulations were performed in the Medical Image Analysis Facility 
at JPL. Most computations were made using a PDP-H/55 Computer System. 

The simulated images were created on a digitized image field of 100 
samples by 100 lines. The distance betv/een picture elements is set at 
1mm of real space. For computational purposes, gray level intensities 
were quantified to 256 discrete levels (8 bits) although the display 
used in this task shows only 16 discrete levels. 

Initial tests were performed on a single point of about 2mm diameter, 
(see Fig. 14A) Line integrals or projections were made at increments of 
7-1/2°. Since the test object was rotationally symmetric, each projection 
should have been identical, but because of the coarse digitization, there 
were noticeable differences which were considered acceptable with regard 
to the assumed resolution limit of 2mm. In order to calculate an enhance- 
ment function later in this process, these projections were averaged 
together to give the following (projection A): 
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Figure 14. Computer Simulated Tests on a Single Point; (a) Simulated 
Object, (b) Reconstructed and (c) Enhanced 
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Without using enhancement, each of the projected lines were then 
used as the source for the reconstructed image in back projection to 
give the image shown in Figure 14b. 

There are some asymmetric characteristics in this reconstruction 
which are artifacts caused by digital roundoff. They can be corrected 
but the increase in numerical manipulation did not appear warranted. 

This image was reprojected and the projections were aga.in averaged 
to give the following (Projection B): 



The Fourier transform Projections A and B were calculated. The 
ratio of the amplitudes of the transforms (at each frequency) was made 
and the resulting ratio was used as input (with zero phase) for an 
inverse Fourier transform back to real space. The results were surprisingly 
brief (Projection C): 
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This curve can bemused as the enhancement kernel for •convolution with each 
projection line prior to reprojection to compensate for the approximate 
|r j ^ weighting which is implicit in the original reconstruction. Figure 
14c shows the effect of applying this correction to the projections used 
to create Figure 14b. As was expected, along with the sharpened central 
peak, the spaicse sampling in the number of projections has caused the 
spikes to appear more prominently. 

A more complex test object was created with points, simulating 1 and 2mm 
separations. (See Figure 15a) Its normal reconstruction and enhanced 
reconstruction are shown in Figures 15b and c respectively . Normal recon- 
struction does not appear to show either 1 or 2mm resolution whereas the 
enhanced image does distinguish the 2mm separations although the noise 
streaks are also prominent. 

The syringe used in Figures 3 and 4, was simulated with three spheres 
at 2,3, and 4mm distance (see Figure 16). The computer program which 
generated the 'image, its projections, enhancement and reprojection -is 
listed in Appendix C. Reconstruction .was first performed using only 90° 
rotation sampling (at 7-1/2° intervals instead of the complete 180° 
sampling). A severe loss in resolution is evident (Figure 17a) and per- 
sists even if enhancement is performed (Figure 17b) although .there is 
some recovery of detail. 

Reconstruction using the complete 180° rotation is shown in Figure 18a. 

f 

A photometric analysis of the image density recorded on Polaroid film 
in Figures 3 to 13 was made with the following results: 
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Figure 15. Computer Simulation Tests on Separate Points: (a) Test Pattern with 1 
and 2 mm Separations (b) Reconstructed, and (c) Enhanced 
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Figure 16. Computer Simulation of Syringe (as shown in Figures 3 and 4) 
Prior to Scan and Reconstruction 
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Figure 17. 90° Potational Scan Beconstruction 



(a) Normal Reconstruction, (b) Enhanced 
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As the film receives more light, (as is used in one mode to represent less 
transmitted sound energy or a denser object), the image becomes less 
responsive to increases in exposure. The result is, that for objects 
which absorb (or reflect) sound energy, there will be a decreased sensi- 
tivity. As less of the input beam is seen by the sensor, then the film 
tends to saturate. Therefore, small differences in a dense object are 
washed out while background striations are visible (See Figures 18b and 4). 
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Enhancement sharpens this test object somewhat (see Figure 19a) and 
simulation of the film characteristics is shown for this enhanced image 
in Figure 19b. Again, the background is enhanced to the detriment of the 
detailed structure of the object of interest. 

The polarity of the camera can be reversed, i.e., increased sound 
density provides increased film exposure (black). Figure 20a is the 
polarity reverse of Figure 18a. Figure 20b now shows a flat background 
and the three spheres are reasonably distinct. 

Figure 21a is the polarity reverse of Figure 19b while Figure 21b is 
enhanced. Although the image is somewhat noiser, the objects are 
delineated more sharply. 
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Figure 20. Reversed Polarity of Figure 18 (a) Normal Reconstruction (b) Polaroid Simulation 
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Figure Z1 


Reversed Polarity of Figure 19 (a) Enhanced Reconstruction (b) Polaroid Simulation 
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Discussion of Computer Results 

The computer has been used to vary parameters and perform 
linear (additive as well as subtractive) and nonlinear operations 
and summations which are difficult to obtain through analog 
manipulation. 

A major parameter is the redundancy of the recorded signal 
which occurs in overlapping scan regions. Also implicit are the 
gaps in information created by the large angular (7-1/2°) incre- 
ments in the scanning sequences. In theory, with an infinitesimal 
angular increment, the redundancy factor approaches lrj"\ In 
order to correct this error it is possible to enhance detail by 
convolving the reciprocal |rj with the resultant final image. 

For the sake of removing noise it turns out that leaving the un- 
enhanced image alone may appear more aesthetic and less distracting, 
but clinically important high frequency detail may be buried and 
lost in such a display. When the enhancement convolution is per- 
formed to recover fine details, not only is random high frequency 
noise made more visible but the error created by undersampling in 
angle is also made visible as “spokes" seen emerging from sharp 
points. 

Enhancement can be performed very rapidly one-dimensional ly 
by convolution of each projection with a one-dimensional pro- 
jection of the |r[ convolution operator. For this work rather 
than use the theoretical factor of jr|, an experimentally 
evaluated enhancement factor was derived. Surprisingly, this 
function was determined to be a brief series of only three 
values -1.2, 3.4, -1.2 as seen in Projection C (p. 23). This 
function, when convolved with each scanned projection integral, 
provides the enhancement required to eliminate the redundancy 
of the low frequencies. 

The results are as expected. It is not clear that 1 mm 
separations are distinguishable while separations of at least 
2 mm are quite clear. (See Fig. 15C - p. 25) 
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Simulation of the nonlinear photometric response of Polaroid 
film also behaved predictably. Image detail whose brightness 
falls in the "toe" or saturated region of the film is lost. 

Digital memory as used in CAT or as simulated here provides 
greater dynamic ranges. 

Finally, the ease of simulation illustrated in this effort and 
the recent availability of inexpensive digital microcomputers 
brings us to consideration of the next step - that of performing 
this entire operation with a real time microprocessor. 
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Additional Considerations 


Time Delay Variations 

In geometric terminology, the path’ or paths which ultrasound 
makes in passing through an inhomogenous medium are geodesic lines. 

It can usually be inferred that if there exists a curve of shortest 
traverse between two points in a medium, then it is geodesic. Geo- 
desics can also be curves of greatest traverse, but in tomogram re- 
construction we are interested in the minimum solution to the varia- 
tional problem of energy passage. 

Straight lines are geodesics for ultrasound only when there is 
no density change capable of bending the sound ray. If one reassembles 
a tomogram from integral projections and assumes the geodesics are 
straight lines when they are not, then the reconstructed image will be 
geometrically distorted. 

JPL has had considerable experience in the reconstruction of trans- 
mission images using the integral projection as a measure of tissue in- 
teraction with transillumlnating ultrasound. Successful transmission 
images have been obtained on in-vitro as well as in-vivo subject matter. 
This experience has shown that if the ultrasound energy corresponding 
to the shortest time of passage is used, and if subsequent signals cor- 
responding to later, less direct paths are removed from the data, then 
reliable images result. 

Consequently, the same type of signal processing was used in this 
tomogram task. The total signal energy density which arrives at the 
receiving transducer is processed and only those components taking a 
minimum time delay in passage, with a time acceptance of one microsecond, 
were used for image reconstruction. In order to use existing software, 
it was necessary to preset the time delay to one particular value for 
all angle positions of the test object relative to the transmitting and 
receiving of ultrasound probes. 
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For transillumination through one direction only, for example, 
imaging through the arm, this preset time delay and window of acceptance 
proves quite useful. It was for that purpose that our existing software 
was designed. However, when reconstructing multiple views taken at all 
possible angles for sound passing through highly elongated objects, this 
existing software is inadequate. 

As an example: If within an irregular mass there are regions of 

differing acoustic velocity, the time delay of the earliest sound will 
vary about any previously preset time delay. If there is an elongated 
region of differing velocity within the mass, then the sound passing 
through the longest axis of that region may arrive at a substantially 
different time than the- sound passing through the shortest axis of that 
region. That is what we observed when testing several i n-vi tro samples 
in which there was a path length variation of greater than 10 centimeters 
over the various angle projections and a tissue velocity considerably 
different than that of the fluid medium in which the specimen was placed. 

In such cases it is impossible to use a fixed-time delay for isolating 
direct from less direct ultrasound, and the tomographic reconstruction 
was not considered useful. 

Refraction of sound did not pose the problem. The problem was due 
to extreme variations in arrival time, such that the energy came in and 
out of the one-microsecond time window. The solution to this difficulty 
was the use of an adaptive time delay which selects the earliest arriving 
sound energy and provides a time window around that early sound. If a 
pulse ultrasound system is used for this method of tomography, the simplest 
circuitry that one can use would be a threshold detector that enables a 
monostable timing circuit for opening and closing the time gate. 

This software modification to optimize display by JPL's specialized 
instrumentation was not attempted during the course of this task since 
to do so would, in the experimenter's judgement, depart from the stated 
objective and divert funds from the goals that had been established. 

It has been adequately demonstrated during the course of this task that 
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detecting the first ultrasound signal and gating out subsequent signals 
provides sufficient data to reconstruct a usable tomogram. This method- 
ology can be applied to existing pulse systems as well as the more power- 
ful phase coherent systems. 

Reconstructing Medium 

Both the analog experiments and computer simulation, to be shortly 
described, have demonstrated that a simple arithmetic process produces 
usable tomograms. The principal difficulties with photographic film as 
a reconstruction medium are the lack of subtractive capabilities of 
photographic emulsion and the inherent nonlinearity of the exposure/ den- 
sity curve. 

The nonl inearl ity of exposure/dens.ity can be overcome by proper gain 
•settings for adj'ustment of contrast. But film can only add light values, 
never subtract them. The result of this is the generation of a background 
bias level on the film which, in order to lie well above the knee of the 
exposure/ density curve, must reduce the usable image contrast. This is 
a problem encountered in X-ray tomography which uses a film medium. 

What is needed is a matrix array which can represent picture elements 
of the reconstructed tomograph and on which the values can be subtracted 
as well as added. When this present task was first contemplated, it was 
not evident that there was going to be a dramatic cost reduction in 
simple microprocessor memory. Recent history has shown that the cost of 
memories may drop sufficiently to allow them to replace film as a simple 
low, cost reconstruction medium. 

The major assertions of this present task have been verified. Namely, 
that a simple mechanical scanner and straightforward arithmetic means can 
be used to generate a less than perfect, but low cost and usable tomogram. 
Furthermore, this arithmetically constructed tomogram compares favorably 
with reconstructions- based on far more extensive algorithms. The sugges- 
tion that a discrete element memory can replace photographic film is no 
retreat from this assertion. 
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APPENDIX A 

RECONSTRUCTION OF THE INTERNAL DENSITY STRUCTURE 
OF AN OBJECT BY ARITHMETIC OPERATION ON A SET OF ITS PROJECTIONS 


There are two well-established computational methods of solving 
the problem of reconstruction from projections. Both of them are 
computationally lengthy. What will be described is a third method. 

A simplified example of the problem is diagrammed in Figure A-1 . 



Figure A-1. Test Problem 

Some object with a distribution of density throughout space, 
f{X^, X 2 ), istrans illuminated by energy E. The transmitting means, 

T, has a radiation of E^, but the absorption by f{x) reduces this 
to E^ by the time it is received by the receiver R. 

The problem is to determine the spatial distribution (f(Xj^, X 2 * * • 
by ^n appropriate measure of the possible projections of energy 
through that distribution. 
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The arithmetic method for accomplishing this is deceptively simple 
and is diagrammed in Figure A-2. 



Figure A-2. Initial Reconstruction Logic 

The energy received through the object is the integral of the density 
along some path, S. Establish a reproduction geometry with a line S 
that is geometrically related to S through the fixed coordinates of 
the density distribution f(Xj, Xg"'). At each place along S add the 
numerical value E.. For every possible transmission path t, establish 
a reproduction line T and add the value of the line integral along t 
as a number to every point on the line T. In the limit as all possible 
views are added, the values on the reproduction plane approach a 
constant plus a close approximation to the distribution that, is desired. 
It is only necessary to add arithmetically positive values until you 
are through with the number of projected views, then the constant 
background level. 
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Solution Using the Radon Transform 


We will be extracting from GEL'FAND (1) for this derivation. 

While the projection operation is not considered in GEL'FAND's 
mathematical treatise, it is apparent that the integral of a 
function over hyperplanes is a generalization of what the scientist 
measures as a projection. 

The functions to be analyzed are infinitely differentiable and 
rapidly decreasing functions on a real affine space. (For simplicity 
the function defined on an n-dimensional space will be noted 

f(X) = f(X^, ^2“'^n^- 

The equation of a hyperplane in this space is given by the inner 
product equality 

(SjX) = Xj + X2 ^ ^ 

where p is a dimensionless constant. 

The volume element is given fay the differential form 

dX = dXi---dX 
1 n 

If a hyperplane of n-1 dimensions is cut through the n-dimensional 
space, the integral of f(X) over the (c,X) = p hyperplane is defined as 
the Radon transform and is given by 

/(c,p) = jf M 6 dX 

where 6 ( ) is the impulse "function" and integration is over the 

whole space. 
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The Radon transform is thus a function defined on a set of 
hyperplanes and is an integral operation on the original function. 

If, as in Figure A-2, the space is two dimensional, the n-1 
dimensional hyperplane is a line S. The integral of f(Xj^*“) 
along that line is a Radon transform of the original function. 

If the original function represents a density, then the integral 
of that density is what is measured by the linear projection. In 
order to provide a physical interpretation of the meaning, consider 
Figure A-3. 



Figure A-3. Projection Correction 

The number p designates the numeric of the hyperplane and, in the 
case of a line, indicates the "parallel location" of that line. This 
is shown in Figure A-3 by the p-axis drawn perpendicular to the pro- 
jection line S. Incrementing p to create the new line S' amounts to 
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a parallel displacement operator. The value of the integral projection 
at every point p is the Radon transform 

The coefficients 

?lXi + " P 

in Figure A-3 define the "angle" of the line with respect to the X-j 5 X 2 
coordinate. These coefficients have the reciprocal dimensionality 
of the coordinates X. 

The dimensionality relations in the inner product representation 
for the Radon transform hyperplane is no accident. The Fourier kernel 
is also defined on an inner product basis and it can be readily shown 
that there is a close tie between the Radon and Fourier transform. 

For example, a Fourier transform in n-dimensions reduces to the Radon 
transform followed by a one-dimensional Fourier transform. In some 
programs this fact can be computationally advantageous since the Radon 
transform has the geometrical meaning of an integral projection. 

In a real affine space, the inverse Radon transform has a peculiar 
characteristic that does not exist if a complex space is used, namely, 
the transform depends upon whether the space is of odd or even dimension. 

The types of problem of immediate interest to the reconstruction 
we are addressing are expressed in two-dimensions, such as Figure A-3. 

t 

If the space is two dimensional the inverse Radon transform 


formula is, 

^ 2 tt /p [g(s)] oiU) 

where. 

CO 

and 

G(5) . ( Ili^ "P 
{p-(c,x)r 
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Other than the statement that the value of the inverse transform 
at some point X depends on the integrals over all possible lines, 
GEL'FAND drops it there. 

Taking a clue from PAPOULIS' (2) discussion of the generalized 
function t"”’, we can proceed as follows. We first find what G(g) 
means in practical terminology. Obviously G(s) involves a generalized 
function, or distribution. 


The generalized derivative of the natural logarithm of the absolute 
value of a coordinate exists and is 


i. ~ d An-lp] 
P dp 


Thus the integral of Injpj with a test function (p) exists and defines 
the distribution Injpj as. 


J ^n|p| <!'(p)dp 


It follows then that the generalized derivative of lnjp| satisfies 


I 


— CO 


«i>(p) dp = 
P 


Continuing, the form for the 



(P) dp 


kernel is' 


CO 


I 


—CO 



dp 


CO 



(m) 


dp 


The value of the function G( ) is then tied in with the integral 
of the second derivative of the Radon transform weighted by a 
logarithmic term. 


Consider the geometry shown in Figure A-4. 
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V 

The signal we record as f is the line integral of the density 
f(X). The p-axis is the axis perpendicular to our illumination path. 
The operator 


_8 

ap 


is the infinitesmal parallel' translation operator for the illumination 
line. 


The value of at the line defined by (?,X) = p^ = {?»Xp} 

is contained in translation of that line as part of 

d y ( 5 ,X) = f (X ) + {other f along the line} 
dp f 

th ^ 

The p second derivative of f then means. 



/ (?.X) = 


dp (Xo) + ^ {other f along the line} 


The term we call G {?) then refers to the rate of increment due to 
f{X^) as we advance the scanning line along the p-axis. This is not 
exact, however, because of the lnjp| weighting term. Now let's take a 
look at that in Figure A-5. 
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At each point in the object, X^, we can consider G (c) to relate 
the inversion of the Radon transform as it relates to rays passing 
through and around but defined by p. In particular the ray passing 
through X^ is the origin of the p-translation and is where (5,X) = p. 
The weighting function becomes ln|p-(c,X)| and has the shape shown in 
Figure A- 5 . 

To first approximation the logarithmic singularity "pulls out" 
the function from inside to outside the integral and "smooths" the 
values immediately around the coordinate thus extracted. This, 
coupled with the face that the log of a big number gets big slowly, 
means we can approximate G (5) by 

^ dp ^ ^^0^ dp C'^'^oo^hed" f(X^) + constant 
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Now the differential (0(5) has an interpretation in ? -space that 
can be inferred from Figure A-6. 



Figure A-6. Description of Coordinate System in ^-Space 

The coefficients, define the "angle" of the scan line relative 
to the coordinate axis of f(x). The meaning of o)(c) is the differen- 
tial for an accumulation along a particular "direction". If we take 
all possible angles, we accumulate along a closed contour, r, which 
encloses the origin in the 5 -space. In a practical sense we only 
need accumulate over one half of r since the projection of a view 
taken at 180° and 0° is the same for most physical problems. 


The significance of 

/(X) - ^ »(5) 

is that of adding up all possible angular views which pass through X. 
Taking the first term in our approximation, 

S-j i 
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we can reconcile the term inside the integral by noting that dp is a 
parallel translation without angle change. In Figure A-6 we can think 
of dp as being proportional to u)(5) since a small angle change from 
one projection to another "adjacent" projection produces a dp change 
of proportional amount. If you v/ill, the dp of the second projection, 
considered itself as a projection on the dp of the first, has its 
orthogonal component (the part not a portion of the first dp) modified 
as the sine of a very small angle change between the two views. 

Thus we can infer that 

In other words, if we take the projection views, and we add them 
up point by point as we look at all possible view angles, the thing 
we accumulate has a functional space dependence that is proportional 
to the original density distribution. 

There are two "error" terms which show in this derivation. One 
is the "smoothed" version of the density and is due to the logarithmic 
kernel. The other is a constant "background". 

The second term is trivial and merely means we must subtract a 
constant value from each image point. 

The first term is more serious and implies we have a slightly out 
of focus image, since the polarity sign is the same for it as the 
proper density. 

In order to estimate this defocussing, we can return to the 
generalized derivative of Injpj 

CO 

dp = - /£n|p|^.^^^p)dp 

CO 00 


I 
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This Is a Hilbert transform relationship defined at the p = 0 
place. The Hilbert transform of a function g (x) is defined as 

'■(X) = V ■’/ 

where P means the Cauchy principal value is to be used (pre-general i zed 
function notation). We can get a reasonable idea of the smoothing 
effect of Injpl by noting how the Hilbert transform of a function 
differs from the derivative of that function. 

Obviously the smoothing effect on an impulse is r”^ where r is 
the radial distance away from the coordinate of the impulse. Discrete 
image points are not impulses, however, they are more like Kronecker 
deltas with finite height. The smoothing effect in that case is very 
small. Other examples, such as edge transitions can be likewise 
exami ned . 

In summary, the method of arithmetic reconstruction appears 
analytically valid and offers a simple solution. 
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Heuristic Argument in Signal Space 

One can offer a reasonable explantion as to why arithmetic 
reconstruction works. Consider Figure A-7. 




Figure A-7. Arithmetic Reconstruction 

The projection, Pj, through some point (X, Y) has a portion of 
its density change, a, due to that point. The integral of the density 
has two terms: a, due to the line integral of everything along the 
path except a, and a itself. 

If lines Lj^, Lg, etc. were laid down on a reproducing plane in 
direct geometric relation to the recording projection lines, and if 
the density of values along the L's were equal to the projection integral, 
the net density at the point A, corresponding to the real a, becomes 

+ £2 (^3 + ^4 +••'+ ^^2n-l ^ ^2n ^ 


50 



77-72 


The total area of the original density, calculated from rays 
drawn from a is 

AREA = + -^2 ^ ^4 +**•+ ■'^2n-l ^ ^2n ^ “ 

In the area computed from a finite number of rays originating 
at a, each ray is added and a is itse’lf only added once. In the 
reconstruction from projection integrals the net value at a is made 
up of the sum of the rays plus N- times the sum of a. Consequently 
the net density at each point (X, Y) is 

D{X,Y) = AREA + (N-1) f(X,Y) + (N-l) (smoothed f(X,Y}) 

where f(X,Y) is the original density distribution. This is exactly 
what we predicted from the inversion of the Radon transform. The 
smoothed term comes about because we assumed a finite projection 
width in Figure A-7. The "star" patterns of converging rays overlap 
very close to (X,Y) and cause a small amount of defocussing. The 
"area" term is the constant background that .must be removed. 

The amount of smoothing is due to the finite resolution of a 
practical system of making projections. It is in fact the amount of 
edge smoothing we will see in any single projection image. If the 
line width of projection becomes vanishingly small, the smoothing 
term will be reduced similarly when a finite number of .projections 
is used. Apparently, the use of an infinite number of projections 
will again lead to a small amount of smoothing. What appears to be 
happening vn'th a finite number of infinitely narrow width projections 
is that the term we call "area" cannot be the true area because there 
are density portions missed between the "spokes". The fluctuation of 
this "area" as a function of (X,Y) is where the smoothing occurs in 
this case. 
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Argument in Fourier Space 

There is yet another way of looking at the method of arithmetic 
reconstruction. The Fourier transform of a projection is the corres- 
ponding slice of the two dimensional Fourier transform of the complete 
object. See Figure A-8. 



P (u, v) = p (x, 


Figure A-8. Reconstruction in Fourier Space 
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If a number of projection views 
two dimensional Fourier transform is 


are taken at different angles, the 
approximated as shown in Figure 9. 



•Figure A-9. Fourier Test Case 

The Fourier transform is not only two dimensional but complex, 
having a magnitude and phase angle as a function of spatial frequency. 
Points, such as (u,v), which lie between the Fourier "spokes" may be 
approximated by extrapolation to those known projection transforms which 
are nearest. As one passes farther out in frequency, such as to point 
(UpVi), the approximations become less accurate per unit increment of 
spatial frequency. 

All Fourier slices pass through the origin. This "overspecification" 
means that the reconstructed image made from a straight summation of 
slices will have a large zero frequency component. 

The lack of inter-spoke preciseness may be interpreted as giving a 
high frequency smoothing. If we simply arithmetically add up all the 
inverse Fourier transforms of this pattern we can expect the resultant 
to be 

F(X,Y) + smoothed f{X,Y) + constant 
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Again, this is exactly what we expect. 
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Appendix B 

Reconstruction in the Presence of Boundary Reflections 


Sound energy will be diminished when passing through a boundary 
separating regions of differing acoustic impedance. Thus, even if 
the regions through which the sound passes have no inherent absorption 
of acoustic energy, there will be a reduction in transilluminated sound. 
This raises the question whether a reconstruction of internal density 
structure from integral projections will be adversely affected by such 
boundaries. 

The presence of such boundary conditions will not, to first order, 
detract from the reconstruction. Consider the idealized case shown 
below. A generalized boundary, B, separates a medium with acoustic 
impedance 2 -^ from a medium with acoustic impedance Z 2 * The boundary 
is assumed to be very thin, since it represents a transition in 
acoustic impedance. 



Fig. 1 
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This- boundary can be considered a continuous locus of scattering 
point centers, much like a densely packed string of very small beads. 

In forming all possible views by sets of projections, each point on the 
boundary will act as a point scatterer independent of all other boundary 
points. In effect the boundary is made up of a contiguous set of point 
scatterers, each of which contributes a constant reduction of passing 
sound energy. The area within the boundary will therefore be reconstructed 
with an indication of a slightly higher average density than would be the 
case had the boundary not been present. 

If linear reconstruction methods are used, then each point scatterer 
will appear in its proper geometric place, have the proper value of 
indicated attenuation, and be out of focus by the same amount due to the 
finiteness of transilluminated beam. This means that the boundary will 
appear in its proper place. 

If the region within the boundary has an additional attenuation 
which represents a different density structure than the region outside 
the boundary, the reconstruction will proceed unmodified by the fact 
that there is a boundary. This is because the sound wave is diminished 
by a constant value, representing the boundary point. 
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Appendix C 

Program Used for Computer Simulation 

FORTRAN IV-PLUS VOl-03 09 22 48 20-SEP-77 

CAT125 FTN /TR BLOCKS 


ORIGINAL PAGE iS 
OP POOR QUALITY 

PAGE i 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 


0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 


DIMENSION CARG<24>! SARG<24) 

INTEGER*2 DTX ( 26. 1 50 ) , ODN. SPAR ( 5 ) 

BYTE 07(100, 100), A ( 5176 ) , LABEL( 5 12) 

DATA LABEL/70*32i "S', 'O', 70*32, "A", "L 368*32/ 

DDN=1+IUNIT(5) 

CALL PARAM(NP, ODN, 1, , " ENTER OUTPUT DDD 
2 FORMATUH , 6512) 

CALL MVL("CAT OUTPUT P ICTURE " , LABEL < 74 ) , 18) 

NT=1 

DO 35 N=l, 100 
DO 35 M=l, 100 
35 DT(N, M)=0 

DO 34 N=24 
DO 34 M=l, 150 
34 DTX(N, M)=0 

DO 400 NR0T=1,24 
ARG=(NR0T~1 )*3 14159/24 
CARG < NROT ) =COS < ARG ) 

400 SARG(NROT)=SIN(ARG) 

DO 1 N=l,100 
DO 1 M=l, 100 
I VW= < N-50 ) **2+ ( M-50 ) **2 

IF((IVW LT 155) AND. ( iUW. GT. 56) )DT(N, M)=a 
I VX= < N-30 ) **2+ ( M-50 > **2 
IFUVX LT 25)DT(N, M)=9 
I VR= < N-69 ) **2+ ( M-50 ) **2 
IFdVR LT 25)DT(N, M)=9 
I VQ= < N~50 ) •*■*2+ ( M-68 ) **2 
IF(IV(5. LT 25)DT(N,M)=9 
1 CONTINUE 

CALL PARAM(NP, NTE3T, 1, , ' TYPE 3 ENHANCE, 20R 3 GEN, OR 1 " 

C GENERATE PROJECTIONS 
DO 500 M~2, 99 
DO 500 M=2, 99 
I=M+1 
J=M-1 
K=N+1 
L=N-1 

500 DT(J,L>= 4*DT(M,N)-r- 1 *( IV(D1 (M, K) )+IV(DT(M, L) )-i-IV(DT( 1, N) ) + 

1 IV(DI ( J, N) ) )+ 05*( IU(DT( I, K) )+IV(DT( I, L) ) + IV(DT( J, K) ) + 

2 IV<DT(J-L))) 

DO 5011 M- 2,99 
DO 5011 N---2, 99 
J=-101-W 
J J“ 1 00" N 
K=iOi-M 
L==^iOO~,M 

5011 DT(K, J)=Dr<L. JJ) 

JF(NTEST EO 1) GOTO 1000 
DO 25 MR0T=1,24 
DO 25 NINC=1, 120 
DO 4 NNINC=1, 120 

Kl-~ ( NN I NC-50 ) *SAR G < NROT ) +50. 5+ ( N I NC-50 ) *C ARG ( NROT ) 
M=<N1NC~50)*SARG(NROT)+50 5 + <NNINC-50)-*CARG(,NR0T ) 

IF(((M GT. 99) OR. (H LT 2>) OR ( (N. GT 99). OR (N. LT 2) ) )G0T04 
D1 X(NROT, NIMC) -DTX (NROT, NINC)+DT(N, M) 
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0054 

0055 

0056 

0057 

0058 

0059 

0060 
0061 
0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 
0081 
0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 
0101 
0102 

0103 

0104 

0105 

0106 
0107 


4 CONTINUE 

25 CONTINUE 

DO 42 N-1, 120 
DTX(25, N>=0 
D041 M=li24 

41 DTX(25, N)=DTX(M, N)+DTX(25, N> 

42 DTX(25, N)=0TX<25, N)/24 
PRINT 30 

DO 5 N =1, 25 
PRINT 28, N 

28 FORMATdH ,12) 

5 PRINT 27, (DTX<N. M), N=l, 120) 

27 FORMAT < IH ,3014) 

C ENHANCE PROJECTION 

IF<NTEST. EQ. 2)<50T017 
DO 200 N==l,24 
DO 201 M=2, 149 
L=M~1 
K=M+1 

201 DTX(26. M)=3 4»-DTX(N, M)“l 2*< DTX < N, L > +DTX < N, K ) ) 

DO 200 M=l, 150 
200 DTX(N, M)=DTX<26, M) 

PRINT 30 

30 FORMAT <1H ) 

DO 501 N=l, 25 
PRINT 28, N 

501 PRINT 27, <DTX(N, M), M=l, 120) 

C REPROJECT BACK 
17 DO 7 N=l, 100 

DO 7 M=l, 100 

7 DT(N, M)=0 

DO 8 N=l, 100 

DO 8 M =1, 100 

DTQ=0 

DO 81 NROT =1,24 

NINC= 50+(N~50 5)*CARG<NRDT) + (M-50. 5)»SARG(NR0T) 
IF<<NINC LT 1) OR <NINC. GT 150))G0T081 
DTQ=DTQ+DTX (NROT, NINC ) 

81 CONTINUE 

8 DT(N, M)=DTQ/100 
PRINT 30 

1000 DO 100 N=l,100 

100 PRINT 2, (DT(N, M), M=20, 84) 

CONTINUE 

CALL OPEN (A, 2560, 1,1, ODN) 

SPAR(1)=100 

SPAR(2)=100 

SPAR(3)=8 

SPAR(4)=4 

SPAR(5)=512 

CALL PLABELC A, SPAR, LABEL) 

DO 150 1=1,100 
CALL PUT (A, 0, lA) 

150 CALL MVL<DT(1> I), A( IA+1), 100) 

CALL CLOSE(A) 

END 
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